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Abstract 

A model of two coupled Ablowitz-Ladik (AL) lattices is introduced. While the system as a whole 
is not integrable, it admits reduction to the integrable AL model for symmetric states. Stability 
and evolution of symmetric solitons are studied in detail analytically (by means of a variational 
approximation) and numerically. It is found that there exists a finite interval of positive values of 
the coupling constant in which the symmetric soliton is stable, provided that its mass is below a 
threshold value. Evolution of the unstable symmetric soliton is further studied by means of direct 
simulations. It is found that the unstable soliton breaks up and decays into radiation, or splits into 
two counter-propagating asymmetric solitons, or evolves into an asymmetric pulse, depending on 
the coupling coefficient and the mass of the initial soliton. 



1 Introduction 



Dynamics of solitary waves (that we will refer to as "solitons" , without implying mathematical rigor) 
in nonlinear lattices is a vast field of research, which is of great interest in its own right and finds 
important physical applications. A majority of nonlinear discrete systems are nonintegrable. How- 
ever, there are paradigm models that are integrablc by means of the inverse scattering transform, 
most notably, the Toda lattice [1] and Ablowitz-Ladik (AL) [2] chain. A nonintegrable general- 
ization of the latter model, in the form of its combination with the discrete nonlinear Schrodinger 
(NLS) equation,was also studied in detail [3, 4]. Although the AL model finds fewer direct appli- 
cations than its discrete-NLS counterpart, it may describe, for instance, ladder-structured lattices 



The availability of exact soliton solutions in integrable models makes it possible to study nontrivial 
dynamical effects in more complex systems that can be built on the basis of integrable ones. In 
particular, a challenging issue is to study solitons and their stability in two coupled systems, which 
are integrable in isolation. A physically important example is a dual-core optical fiber (the so-called 
directional coupler [6]). While each core is described by the integrable NLS equation, the coupled 
system is not integrablc, and it gives rise to a new effect: an obvious symmetric soliton solution, 
with its energy equally split between the cores, becomes unstable if the energy exceeds a certain 
critical value [6]. As a result of the onset of the instability, a pair of new solitons with a broken 
symmetry emerge (they are mirror images to each other). Later, similar bifurcations destabilizing 
symmetric solitons and replacing them by asymmetric ones were found in a dual-core system with 
quadratic (rather than cubic) nonlinearity [7], and in a dual-core fiber with the Bragg grating [8]. 

A natural step is to introduce a system of two coupled AL chains and investigate the stability and 
nonlinear evolution of its solitons. A system of two AL chains with a coupling that does not admit 
reduction to the usual AL model was introduced in Ref. [9], and solitary waves, as well as moving 
breathers, were found in it. In this work, we focus on a system of coupled AL chains that admits 
reduction to the integrable AL model and seems, as a matter of fact, more natural. The system is 



where w„ and Vn are complex dynamical variables at the n-th site of the chain, the overdot stands 
for the time derivative, and e is a real coupling constant. Under the symmetric reduction Un = Vn, 
Eqs. (1.1) and (1.2) reduce to the AL model proper. 

Results presented in this work demonstrate that the system of equations (1.1) and (1.2) is definitely 
nonintegrable, unlike the AL model. The present model conserves only two dynamical invariants, 
viz., the Hamiltonian, 



[5]. 



iUn + (1 + \Unf) [{Un+1 + n„_i) -|- e(v„+i -I- ^^n-l)] = 0, 

iVn + (1 + l^nP) [{Vn+1 + Vn-l) + €{Un+l + ti„_l)] = 0, 



(1.1) 
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and the total "mass" , 
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(note that the masses corresponding to the u- and i;-fields are not conserved separately). 

The symmetric reduction, Un{t) = Vn{t), generates a solution to Eqs. (1.1) and (1.2) which is 
tantamount to the well-known AL soliton, 

Un{t) = Vn{t) = (sinha) sech (an + vt — Oq) exp [—i {hn + ujt — (f)Q)\ (1.5) 

(in this paper, we only consider bright solitons). Here, 9q and are position and phase constants, 
while parameters a and h determine the soliton's velocity and frequency, 

f = 2 (1 + e) sinha • sin 6, a; = — 2 (1 + e) cosh a • cos 6 (1-6) 

[notice that the coupling parameter e does not appear in the expression (1.5) as it is absorbed into 
the time variable, which amounts to a rescaling of the velocity and frequency in Eq. (1-6)]. The 
value of the mass (1.4) corresponding to the soliton (1.5) is 

Msoi = 4a. (1.7) 

In this work, we primarily focus on the stability of symmetric solitons (1.5) with zero velocity, i.e., 
6 = 00 = 0. 

The underlying equations (1.1) and (1.2) also allow the anti-symmetric reduction, Un{t) = —Vn{t). 
But there is no necessity to treat this case separately, as it is tantamount to the symmetric reduction 
with the change e — e. We stress that, in this work, we consider both positive and negative values 
of e. In particular, the symmetric-soliton solution given by Eqs. (1.5) and (1.6) remains valid also 
in the case 1 -|- e < 0. 

The rest of the paper is organized as follows. In section 2, we put forward an analytical approach 
to the stability problem, based on a variational approximation (a review of the application of this 
technique to solitons was recently given in Ref. [10]). This approach, which we employ in its 
simplest form, will produce partial information on the stability, as it will be seen from comparison 
with results of numerical computation of stability eigenvalues. The numerical eigenvalues are 
presented in section 3. Finally, direct numerical simulations of the full system of Eqs. (1.1) and 
(1.2), showing nonlinear development of the instability, are displayed in section 4. 



2 The variational approximation 

A feasible source of instability of the two-component soliton (1.5) is an eigenmode of small pertur- 
bations splitting the two components. To accommodate this mode, we adopt the following ansatz 
for a perturbed soliton, 

Un{t) = (sinha) sech [a (n — a:(t))] exp [+ic(t)n + i(/)(t)] , (2-1) 
Vn{t) = (sinha) sech [a (n + x(t))] exp [—ic(t)n + i0(t)] , (2-2) 

where 2x{t) is a small time-dependent separation between centers of the u- and f-components, and 
2c{t) is a dynamically conjugate variable, viz., a wavenumber difference between the components, 
while (l){t) is a common time-dependent phase of both components. Note that this ansatz does not 
include another possible perturbation mode, which may introduce a phase difference between the 
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two components of the soliton, its conjugate variable being the ampUtude difference between the 
components [11]. We focus on the restricted ansatz (2.1), (2.2) which accounts for the position 
sphtting, as its counterpart which also takes into consideration the phase differences turns out to 
be cumbersome to calculate, therefore it will not be pursued in this paper. 

To apply the variational technique, we need a Lagrangian of the coupled system (1.1), (1-2), which 
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(2.3) 



where H is the Hamiltonian defined in Eq. (1.3). An effective Lagrangian can be calculated by 
substituting the ansatz (2.1), (2.2) into the Lagrangian (2.3) and expanding it in powers of the small 
separation parameters x and c up to quadratic terms, which are necessary to generate perturbed 
equations of motion linear in x and c. Due to the obvious symmetry of the Lagrangian (2.3), the 
time derivative of x will not appear in the result, while the time derivative c of c may only appear 
linearly (being multiplied by x). Finally, making use of some formulas for infinite sums related 
to the soliton's waveform [involving sech (an)], which were borrowed from Ref. [4]), the effective 
Lagrangian is found in the following form 



where A is a positive constant defined as 

1 



acx + (sinh a) {c^ + Aa^e j 



(2.4) 



A = sinh a ^ < -sech (an) sech (a (ra — 1)) sech^ (an) + sech^ (a (n — 1)) 



+sech^ (an) sech^ (a (n — 1)) sinh (an) sinh (a (n — 1))| — 1 . 



(2.5) 



In particular, ^4 = 1 in the continuum limit (a — ^ 0), and A 6exp (—2a) in the ultradiscrete limit 
(a —>■ oo). 



The effective Lagrangian (2.4) immediately gives rise to an evolution equation 

X -\- ^sinh^ a}j Ax = Q , 



(2.6) 



which predicts that the symmetric soliton (1.5) is unstable against the splitting perturbation mode 
if e < 0. In this case, the instability growth rate is 



A = (2 sinh a) \/—Ae , 
as predicted by Eq. (2.6), i.e., by the variational approximation. 



(2.7) 



3 Numerical analysis of the instability of symmetric solitons 



A general stability analysis of the symmetric soliton (1.5) with v = 9q = Q assumes that a perturbed 
solution is taken as 

Un = {Un + Un)e'-^\ Vn = {Un + t*„)e"''^*, Un = (sinha) sech(an), (3.1) 
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where the frequency u is the same as defined in Eqs. (1-6) (with h = 0), and Un,Vn are in- 
finitesimal perturbations. Numerically simulating equations (1.1) and (1.2) linearized about the 
unperturbed soliton solution, we have found that the most unstable perturbation modes are always 
anti- symmetric ones, i.e., with Un = —Vn (which is not surprising, as the same is true for a mode 
destabilizing symmetric solitons in all the previously studied continuum models of the dual-core 
type [6, 7, 8]). Thus, we look for eigenmodes of the form 

Un = -Vn = fn exp (Xt) + exp {\*t) , (3.2) 

where A is a (generally complex) instability growth rate. Substituting Eqs. (3.1) and (3.2) into 
Eqs. (1.1) and (1.2) and linearizing these equations, one eventually arrives at an eigenvalue problem 
based on the equations 

{i\ + a;) /„ + (1 - e)(l + \Unf){fn+l + fn-l) + (1 + e)Un{Un+l + Un-l){fn + 9n) = 0, (3.3) 

{-iX + u)gn + {l- e)(l + \Un?){gn+i + gn-i) + (1 + e)C/„(C/„+i + Un-l){gn + fn) = 0, (3.4) 

which are supplemented by the boundary conditions demanding that, for discrete eigenmodes, the 
fields fn and gn vanish at |n| = oo. 

The eigenvalue problem (3.3) and (3.4) was solved numerically by the shooting method. Figure 
1 presents unstable eigenvalues as a function of the coupling constant e for three different fixed 
values of the soliton parameter, a = 0.8, 0.9, and 1 [see Eqs. (1.5) and (1.6)]. We have checked 
by simulating evolution governed by the linearized equations that the unstable eigenvalues which 
are presented in Fig. 1 are the most unstable ones. We suspect that these are all the unstable 
eigenvalues in the above system (3.3) and (3.4). Even if they are not, other eigenvalues are less 
unstable and thus less significant. 

Our numerical results for the eigenvalue problem show that, for each value of a, the unstable 
eigenmodes at e < are odd, i.e., /_„ = —fn and g-n = —gn- When e > 0, the eigenmodes are 
even, i.e., /_„ = fn and g-n = gn- Furthermore, when negative e has a small absolute value, the 
inspection of the unstable eigenmode shows that it corresponds to a splitting of the two components 
(appearance of separation between their centers) in the soliton (1.5). Recall that precisely this type 
of the perturbation mode was assumed in the ansatz (2.1) and (2.2). On the other hand, when 
e > and small, the unstable eigenmode corresponds to a phase-difference instability, which is not 
accounted for by the ansatz. 

As is seen in Fig. 1, in the cases a = 0.9 and 1.0 the solitons (1.5) are unstable for all values of e. The 
instability may be both non-oscillatory or oscillatory (corresponding to a real or complex eigenvalue 
A, respectively) in different intervals of e. However, a stability window, 0.605 < e < 1.654, is found 
at a = 0.8. More detailed computations show that the stability window opens up at a critical value 
of the soliton parameter, Ocr ~ 0.881, and persists in the region a < acr- At the point a = Ocr, the 
stability window appears (up to the accuracy of the numerical data) at the value of the coupling 
constant Ccr = 1. 

Inspection of Fig. 1 leads to the following general conclusions. First, all symmetric solitons are 
unstable when e < 0. In particular, near e = — 1, the real part of the unstable eigenvalue is very 
small but positive. Second, narrower solitons (with a larger mass, i.e., larger value of a) are more 
unstable. Third, when the soliton mass is below the threshold value 4acr, the soliton is stable inside 
a certain positive-e window. These features are somewhat similar to those in the above-mentioned 
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continuum models of the nonlinear-optical dual-core systems with the linear coupling between the 
cores [6, 7, 8]. In those models too, symmetric solitons become unstable if their mass exceeds a 
critical value, while the coupling constant is positive, and they axe never stable if the coupling 
constant is negative. However, the stability window observed in Fig. 1 when a < Ocr is different 
from stability regions in the continuum dual-core models, as the latter stability regions have no 
right cndpoint (they are semi- infinite) . 

It is natural to compare the numerically found unstable eigenvalues with the one (2.7) that was 
predicted by the variational approximation in the previous section. This comparison is meaningful 
only when e < and |e| <C 1, as the actual unstable mode corresponds to the variational ansatz 
(2.1) and (2.2) only in this case. For a = 0.8, the comparison is displayed in Fig. 2. It is not 
surprising that the deviation between the analytical and mimcrical results is large for large values 
of |e|, when the coupling term can make the perturbed soliton strongly different from the ansatz 
(2.1), (2.2). For small |e|, the analytical value is close to the numerical one; the remaining difference 
might be due to the fact that, if the two components of the soliton are separated by some distance, 
the symmetry of each component relative to its center is broken by the coupling term, which is not 
taken into regard by the ansatz. 



4 Direct simulations of the instability development 

In order to understand the evolution of unstable symmetric solitons, direct numerical simulations of 
the full nonlinear equations (1.1) and (1.2) were performed. To this end, the initial configuration was 
taken in the form corresponding to the ansatz (1.1), (1-2), with an extra perturbation parameter, 
namely, a phase difference 2V' between the two components of the soliton: 



Different typical outcomes of the instability development can be illustrated by a set of contour-plot 
pictTircs pertaining to several characteristic values of e, while all the other parameters arc fixed. In 
Figs. 3 through 7, we display the pictures for a representative case a = 0.8, xq = cq = ipQ = 0.01. In 
fact, exact initial values of the small perturbations are not important, while the value of the soliton 
parameter a is a significant one. In each figure, contours in the left and right panels represent the 
evolution of and respectively. All the contours start at the level 0.05 and increase with 
an increment of 0.2. 

First of all. Fig. 3 shows that if e belongs to the stability interval, see the upper panel in Fig. 1 
(in Fig. 3, e = 1), the initial perturbation indeed does not trigger any instability. Next, Fig. 4 
shows that, in the case e = — 7, which corresponds to strong instability, the soliton gets completely 
destroyed, decaying into radiation. Decreasing the absolute value of negative e, i.e., proceeding to 
weaker instability, according to Fig. 1, we observed a trend of splitting of the original unstable 
symmetric soliton into two moving ones, each being a two-component pulse. An example of that 
is shown in Fig. 5 for e = —0.08. It is noteworthy that the symmetry of each secondary soliton 
is strongly broken (the amplitude of one component is definitely larger than that of the other), 
but they are (at least, approximately) mirror images to each other, so that the global symmetry is 
conserved. 



Un{t) 
Vn{t) 



(sinha) sech (a (n — xq)) exp {—icQin + iijjQ) . 
(sinha) sech (a (n -|- xq)) exp {+icon — iVo) • 



(4.1) 
(4.2) 



6 



In the case when the symmetric sohton is unstable at positive e, the character of the instabihty is 
quite different from that described above for e < 0. In this case, the instability is non-oscillatory. In 
Fig. 6, which pertains to e = 0.3, a noteworthy result is the formation of a soliton with a strongly 
broken symmetry and (nearly) periodic internal oscillations, which is accompanied by emission of 
small amounts of radiation. This instability-induced spontaneous symmetry breaking resembles 
what is known in the above-mentioned continuum models of the dual-core type [6, 7, 8]. Of course, 
in all the cases when spontaneous symmetry breaking is observed, its sign (e.g., the difference 
between the left and right panels in Fig. 6) is determined by a random initial perturbation. 

Finally, Fig. 7 shows the situation in the case e = 2, which lies beyond the right border of the 
stability interval in Fig. 1 (recall that no such border occurs in the contimium dual-core models). 
As is seen, in this case the instability also breaks the symmetry of the soliton (in the beginning). 
However, an unusual feature here is that the centers of both components oscillate, essentially, in- 
phase. Unlike what was seen in all the other figures, spontaneous onset of such in-phase oscillations 
seems to violate the momentum conservation; however, one should keep in mind that lattice systems 
conserve no momentum, in view of the lack of the continuous translational invariance in them. 
Nevertheless, the appearance of such a dynamical state is a remarkable fact which may deserve 
further investigation. 



5 Conclusion 



In this work, we have introduced a model of two coupled Ablowitz-Ladik chains. While the system 
as a whole is not integrable, it admits reduction to the integrable AL model for symmetric states. We 
have studied the stability and nonlinear evolution of stationary symmetric solitons in detail. Both 
the analytical consideration, based on the variational approximation, and numerical computation of 
the instability eigenvalues have demonstrated that the soliton may be unstable. Numerical results 
also show that, provided that the soliton's mass is below a critical value, there exists a finite interval 
of positive values of the coupling constant e in which the symmetric soliton is stable. Comparison of 
the approximate analytical and exact numerical unstable eigenvalues shows that the agreement is 
reasonable for small negative values of the coupling constant. Evolution of the unstable symmetric 
soliton was further studied by means of direct simulations. It was found that the unstable soliton 
can decay into radiation, or split in two counter-propagating asymmetric solitons, or evolve into an 
asymmetric pulse, depending on the values of the coupling coefficient e and the soliton mass. 
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Figure 1: The numerically found unstable eigenvalue A versus the coupling coefficient e for the 
symmetric stationary soliton at different fixed values of the soliton parameter a [see Eqs. (1.5) 
and (1.6)]. The solid and dashed curves show, respectively, the real and imaginary parts of A. 
The imaginary part is not shown inside the stability window in the upper panel, where A is pure 
imaginary. In the interval around the point e = — 1, the real part of A is very small but finite in all 
three panels. 
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Figure 2: The instability growth rate predicted by the variational approximation, see Eq. (2.7) 
(the dashed curve), vs. its numerically computed counterpart (the solid curve). 
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Figure 3: Evolution of the initially perturbed symmetric soliton in the case a = 0.8, e = 1, when 
the soliton is stable. In this figure and below, all contours start at the level 0.05 and increase with 
an increment of 0.2. 
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Figure 5: Splitting of an unstable symmetric soliton in the case a = 0.8, e = —0.08. 
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